Hyperspectral imaging

Step One – Reading Hdf5 data into R

Hyperspectral data consists of many bands that cover the electromagnetic spectrum.

We will use the raster and the rhdf5 libraries to read in the hdf5 data that contains hyperspectral data for the San Joaquin NEON field site in Domain 17. We also need the h5metadata function which is a custom script that allows us to import metadata including the coordinate reference system into R. We will use this information to plot things and to perform the math.

library(raster)

## Loading required package: sp

library(rhdf5)

f <- "/Users/THART/scratch/ESAData_2014/Part2_Spectrometer/HDF5File/SJER_140123_chip.h5"
h5ls(f,all=T)

##   group        name         ltype corder_valid corder cset       otype
## 0     / Reflectance H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
## 1     /        fwhm H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
## 2     /    map info H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
## 3     / spatialInfo H5L_TYPE_HARD        FALSE      0    0   H5I_GROUP
## 4     /  wavelength H5L_TYPE_HARD        FALSE      0    0 H5I_DATASET
##   num_attrs  dclass          dtype  stype rank             dim
## 0         5 INTEGER  H5T_STD_I16LE SIMPLE    3 477 x 502 x 426
## 1         2   FLOAT H5T_IEEE_F32LE SIMPLE    2         426 x 1
## 2         1  STRING     HST_STRING SIMPLE    1               1
## 3        11                                  0                
## 4         2   FLOAT H5T_IEEE_F32LE SIMPLE    2         426 x 1
##            maxdim
## 0 477 x 502 x 426
## 1         426 x 1
## 2               1
## 3                
## 4         426 x 1

#make sure the hdf5 metadata function is loaded
source("/Users/THART/scratch/neonESA2014/lessons/support/h5metadata.R")

Let’s use our custom function to grab metadata from the hdf5 document.

spinfo <- h5metadata(f,"spatialInfo",11)

Next, let’s read in the wavelength center associated with each band in the hdf5 file. What wavelength is band 19 associated with? (hint: look at the wavelengths vector that we just imported and check out the data located at index 19)

#read in the wavelength information from the Hdf5 file
wavelengths<- h5read(f,"wavelength",index=list(1:426,1))

Now - let’s have a look at one of the bands. Each band in the hdf5 file represents a particular Let’s check out the green band… which is band 34. For fun - ask R to tell you what wavelength value band 19 is.Slice up bands make a quick plot of the data.

b34<- h5read(f,"Reflectance",index=list(1:477,1:502,34))
## Convert from array to matrix
b34 <- b34[,,1]

#let's plot one band, to see what it looks like ... 

image(b34)

#hmmm it look a bit fishy. let's look at the distribution of values
#let's have a look at the range of values in our data.
hist(b34,breaks=40,col="darkmagenta")

hist(b34,breaks=40,col="darkmagenta",xlim = c(0, 5000))

#what's going on towards the right of the plot? It looks like we have some pixels that are skewing how this renders
hist(b34, breaks=40,col="darkmagenta",xlim = c(5000, 15000),ylim=c(0,100))

#from the metadata we know that 15000 is "no data" - let's set it to NA so R doesn't try to render those pixels
b34[b34 > 14999] <- NA

#let's run a simple log on the values to try to factor out those high numbers just a bit. we will learn how to sharpen our data later.
image(log((b34)))

so… now we have an image - but perhaps not quite oriented in the correct direction. This is because R reads in matrices starting from the upper left hand corner. Whereas, most rasters read pixels starting from the lower left hand corner. We can deal with this issue but creating a proper raster in R that will read in pixels as other software like QGIS and ENVI do.

let’s try again but this time, let’s flip the order of things

b34_2<- h5read(f,"Reflectance",index=list(1:477,502:1,34))
## Convert from array to matrix
b34_2 <- t(b34_2[,,1])
image(log(t(b34_2)))

#sweet! that looks a bit more like the image that we want to see

Ok, so now let’s create a meaningful raster, that is plotted in geographic space from our data. to do this we need to know the coordinate reference system, and the location of the location of the first pixel (located in the lower left hand corner of the raster) we also need the resolution or size of each pixel in the data. Now NOTE that we are using the original matrix - b34 that we created above. if you remember correctly, this matrix is flipped along the Y axis. this is OK because the R raster function is smart enough to know to how to read in the data properly (starting at the lower left hand corner!).

b34r <- raster((b34))

## Now we need to grab the extent in a bounding box.  Luckily this is in the metadata,  We'll need it in the form of LL lon, UR lon, LL lat, UR lat

ex <- sort(unlist(spinfo[2:5]))
e <- extent(ex)
extent(b34r) <- e
plot(b34r)

Awesome! Now we’ve plotted one band in geographic space. Try plotting some other bands and see what it looks like. Notice that the images that we plot lack contrast… we’ll get to that next!

Now, let’s try and generate a full color image with the RGB bands. But before we start, let’s write a function to handle some of the basic cleaning we did earlier, that way we can bulk process bands.

# f: the hdf file
# band: the band you want to grab
# returns: a cleaned up HDF5 reflectance file
getBandMat <- function(f, band){
  out<- h5read(f,"Reflectance",index=list(1:477,1:502,band))
  ## Convert from array to matrix
  out <- t(out[,,1])
  out[out > 14999] <- NA
  return(out)
}
band2rast <- function(f,band){

out <-  raster(getBandMat(f,band),crs="+zone=11N +ellps=WGS84 +datum=WGS84 +proj=longlat")

  ex <- sort(unlist(spinfo[2:5]))
# If you want to stay in UTM's you can use the code below and comment out the two lines above. Note that these were calculated externally and not embedded in the metadata.  In the future they will be embedded in the metadata of NEON HDF5 files.

#out <-  raster(getBandMat(f,band),crs="+zone=11 +units=m +ellps=WGS84 +datum=WGS84 +proj=utm")
# ex <- c(256521.0,256998.0,4112069.0,4112571.0)
  e <- extent(ex)
  extent(out) <- e
  return(out)
}



stackList <- function(rastList){
  ## Creates a stack of rasters from a list of rasters
  masterRaster <- stack(rastList[[1]])
  for(i in 2:length(rastList)){
    masterRaster<-  addLayer(masterRaster,rastList[[i]])
  }
  return(masterRaster)
}

rgb <- list(58,34,19)
rgb_rast <- lapply(rgb,band2rast, f = f)
## Add the names of the bands so we can easily keep track of the bands in the list
names(rgb_rast) <- as.character(rgb)
### Check with a plot
plot(rgb_rast[[1]])

rgb_stack <- stackList(rgb_rast)
plot(rgb_stack)

plotRGB(rgb_stack,r=1,g=2,b=3, scale=300, stretch = "Lin")

writeRaster(rgb_stack,file="test6.tif",overwrite=TRUE)

## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
## Path to GDAL shared files: /usr/local/Cellar/gdal/1.11.0/share/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)

## class       : RasterBrick 
## dimensions  : 502, 477, 239454, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 1.159e-05, 8.757e-06  (x, y)
## extent      : -119.7, -119.7, 37.12, 37.13  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/THART/scratch/neonESA2014/lessons/test6.tif 
## names       : test6.1, test6.2, test6.3 
## min values  :     123,     116,      84 
## max values  :   14343,   14896,   13805

a note about image stretching notice that hte scale is set to 300 on the RGB image that we plotted above. We can adjust this number and notice that the image gets darker - or lighter

If you want to play around a bit with this – try plotting the RGB image using different bands. Here are some suggestions. Color Infrared / False Color: rgb: (90,34,19) SWIR, NIR,Red Band – rgb (152,90,58)

More on Band Combinations: http://gdsc.nlr.nl/gdsc/en/information/earth_observation/band_combinations

We can also plot our image on a map of the US

library(maps)
map(database="state",region="california")
points(spinfo$LL_lat~spinfo$LL_lon,pch = 15)

### Add raster

Now the fun stuff! Let’s create NDVI or Normalized Difference Vegetation Index NDVI is simply a ration between the Near infrared portion of the spectrum and the red. Please keep in mind the there are different ways to aggregate bands when using hyperspectral data. We are just showing you how the do the math. This is not necessarily the best way to calculate NDVI using hyperspectral data!

ndvi_bands <- c(58,90)

ndvi_rast <- lapply(ndvi_bands,band2rast, f = f)
ndvi_stack <- stackList(ndvi_rast)
NDVI <- function(x) {
  (x[,2]-x[,1])/(x[,2]+x[,1])
}
ndvi_calc <- calc(ndvi_stack,NDVI)
plot(ndvi_calc)